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Abstract. This paper presents a complete methodology for the Bayesian inference of a 
semi-Markov process, from the elicitation of the prior distribution, to the computation 
of posterior summaries, including a guidance for its JAGS implementation. The inter- 
occurrence times (conditional on the transition between two given states) are assumed to 
be Weibull-distributed. We examine the elicitation of the joint prior density of the shape 
and scale parameters of the Weibull distributions, deriving in a natural way a specific 
class of priors, along with a method for the determination of hyperparameters based on 
"historical data" and moment existence conditions. This framework is applied to data of 
earthquakes of three types of severity (low, medium and high size) occurred in the central 



Northern Apennines in Italy and collected by the CPTI04 ( 2004 ) catalogue. 



1. Introduction 



Markov renewal processes or their semi-Markov representation have been considered in 
the seismological literature since they are models which allow the distribution of the inter- 
occurrence times between earthquakes to depend on the last and the next earthquake and 
to be not necessarily exponential. The time predictable and the slip predictable models 



studied in Grandori Guagenti and Molina ( 1986 ) , in Grandori Guagenti et al. ( 1988 ) and in 



Betro et al. (1989 ) are special cases of Markov renewal processes. These models are capable 



of interpreting the predictable behavior of strong earthquakes in some seismogenetic area. 
In these processes the magnitude is a deterministic function of the inter-occurrence time. 
A stationary Markov renewal process with Weibull inter-occurrence times has been studied 



from a classical statistical point of view in Alvarez (2005). The Weibull model allows to 
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consider monotonic hazard rates; it contains the exponential model as special case which 
gives a Markov Poisson point process. In that paper the parameters of the model have 
been fitted to the large earthquakes in North Anatolian Fault Zone through maximum 
likelihood and the Markov Poisson point process assumption is tested. In order to capture 



a non monotonic behavior in the hazard, in Garavaglia and Pavani (2009) the model 



of Alvarez (2005) is modified and a Markov renewal process with inter-occurrence times 
that are mixtures of an exponential and a Weibull distribution is fitted to Turkish data. 



In Masala (2012) a parametric semi-Markov model with generalized Weibull distribution 
for the inter-occurrence times is considered and fitted to Italian earthquakes. This semi- 



Markov model with generalized Weibull distributed times has been first used in Foucher 



et al. (2009) to study the evolution of HIV infected patients. Votsi et al. (2012) consider 



a semi-Markov model for the seismic hazard assessment in the Northern Aegean sea and 
they estimate the quantities of interest (semi-Markov kernel, Markov renewal functions, 
etc.) through a nonparametric method. 

While a wide literature concerning classical inference for Markov renewal models for 
earthquake forecasting exists, to our knowledge in this context a Bayesian approach is lim- 



ited. Patwardhan et al. (1980) consider a semi-Markov model with log-normal distributed 
discrete inter-occurrence times that apply to the large earthquakes in the circum-Pacific 
belt. They stress the fact that using Bayesian techniques is relevant when prior knowl- 



edge is available and is fruitful even if the sample size is small. Marin et al. (2005) have 
also employed semi-Markov models in the Bayesian framework but applied to a completely 
different area: sow farm management. They use WinBugs to perform computations (but 
without giving details) and elicit their prior distributions on parameters from knowledge 
on farming practices. 

From a probabilistic viewpoint, a Bayesian statistical treatment of a semi-Markov pro- 
cess amounts to model the data as a mixture of Semi-Markov processes, where the mixing 
measure is supported on the parameters, by means of their prior laws. A complete char- 



acterization of such a mixture is given in Epifani et al. (2002). 
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In this paper we develop a parametric Bayesian analysis for a Markov renewal process 
modeling earthquakes in an Italian seismic region. The magnitudes are classified into three 
categories according to their severity and these categories represent the state visited by the 



process. Following Alvarez (2005), the inter-occurrence times are assumed to be Weibull 
random variables. The prior distribution of the parameters of the model is elicited using an 
"historical dataset". The time of the last occurrence in the historical dataset constitutes 
the origin for the "current sample" , which is formed by the sequences of earthquakes and 
the corresponding inter-occurrence times collected up to some time T to form the likelihood 
function. When T does not coincide with an earthquake the last observed inter-occurrence 
time is censored. The posterior distribution of the parameters is obtained through Gibbs 
sampling and the following summaries are estimated: the transition probabilities, the 
shape and scale parameters of the Weibull waiting times for each transition and the so- 
called cross state-probabilities (CSPs). The transition probabilities indicate whether the 
strength of the next earthquake is in some way dependent on the strength of the last one; 
the shape parameters of the waiting time indicate whether the hazard rate between two 
earthquakes of given magnitude classes is decreasing or increasing; the CSPs indicate what 
is the probability that the next earthquake occurs at or before a given time and is of a given 
magnitude, conditionally on the waiting time from last earthquake and on its magnitude. 
The paper is organized as follows. In Section [2] we illustrate the dataset. It is one of 



the sequences analyzed by Rotondi (2010). Section |3| introduces the parametric Markov 
renewal model. Three magnitude classes are considered and a Weibull distribution is 
assumed for the inter-occurrence times. Section [4] deals with the elicitation of the prior, 
which is based on a generalized inverse Gamma distribution, considering also the case of 
scarce prior information. Section [5] contains the data analysis with the estimation of the 
above-mentioned summaries. The matter of whether the observed earthquakes reflect a 
time predictable or a slip predictable model is also discussed. Section [6] is devoted to some 
concluding remarks. An appendix contains the detailed derivation of the full conditional 
distributions and the JAGS implementation of the Gibbs sampler. 
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Figure 1. Map of Italy with dots indicating earthquakes with magnitude 
M w > 4.5 belonging to macroregion MR 3 (Rotondi (2010)). Inclusion in the 
macroregion was based on the association between events and seismogenetic 
sources; the region contour has only an aesthetic function. 



2. A TEST DATASET 



We test our method on a sequence of seismic events chosen among those examined 



by Rotondi (2010), which was given us by the author. The sequence collects events 



that occurred in a tectonically homogeneous macroregion, identified as MR3 by |Rotondi" 



(2010) and corresponding to the central Northern Apennines in Italy. The subdivision 



of Italy into eight (tectonically homogeneous) seismic macroregions can be found in the 



DISS (2007) and the data are collected by the CPTI04 (2004) catalogue. Considering 



earthquakes with magnitude M w > 4.5, the sequence is complete from year 1838: using 

a lower magnitude would make the completeness of the series questionable, especially in 

its earlier part. The map of these earthquakes marked by dots appears in Figure [T] As a 

lower threshold for the class of strong earhquakes we choose M w > 5.3, as suggested by 
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Rotondi (2010). Then a magnitude state space with three states is obtained by indexing an 



earthquake by 1, 2 or 3 if its magnitude belongs to intervals [4.5, 4.9), [4.9, 5.3), [5.3, +oo) 



respectively. Magnitude 4.9 is just the midpoint between 4.5 and 5.3. Rotondi (2010) 
considers a nonparametric Bayesian model for the distribution of the inter-occurrence time 
between strong earthquakes (i.e. M w > 5.3), after a preliminary data analysis which rules 
out Weibull, Gamma, log-normal distributions among other frequently used holding time 
distributions. On the other hand, with a Markov renewal model, the sequence of all the 
inter-occurrence times is subdivided into shorter ones according to the magnitudes, so that 
we think that a parametric distribution is a viable option. In particular, we focussed on 
the macroregion MR 3 because the Weibull distribution seems to fit the inter-occurrence 
times better than in other macroregions. This fact is based on qq-plots. The qq-plots 
for MR3 are shown in Figure [2] The plot for transitions from 1 to 3 shows a sample 
quantile that is considerably larger than expected. The outlying point corresponds to a 
long inter-occurrence time of about 9 years, between 1987 and 1996, while 99 percent of 
inter-occurence times are below 5 years. Obviously, the classification into macroregions 
influences the way the earthquake sequence is subdivided. 



3. Markov Renewal Model 

Let us observe over a period of time [0, T] a process in which different events occur, with 
random inter-occurence times. Let us suppose that the possible states of the process are 
the points of the finite set E — {1, . . . , s} and suppose that the process starts in state jo- 
Let us denote by r the number of times the process changes state in the time interval [0, T] 
and let Si denote the time of the i-th change of state. Hence, < si < • • • < s T < T. Let 
jo, ji, ■ ■ ■ , j r be the sequence of states visited by the process and Xi the sojourn time in 
state jj_i, for i — 1, . . . r. Then 



(1) 



Xi = Si - Sj_i forz = l,...,r 

5 









-2.0 -1.0 0.0 0.5 
theor. quant. 





-2.0 -1.0 0.0 
theor. quant. 




-1.5 -1.0 -0.5 0.0 0.5 
theor. quant. 



FIGURE 2. Weibull qq-plots of earthquake inter-occurrence times (central 
Northern Appennines) classified by transition between magnitude classes. 



with s : = 0. Furthermore, let ut be the time passed in j T i.e. 



(2) 



ut = T — s T , 



Time ut is a right-censored time. Finally, our data are collected in the vector (j,x, ut), 
where (j,x) = (j„,x n ) n=1> ... ;T . 

In what follows, we assume that the data are the result of the observation, on the time 
interval [0, T], of an homogeneous Markov renewal process (J n ,X„) n > starting from j . 
This means that the sequence (J n ,X n )„> satisfies 



(3) 



P(J =Jo) = l, P(X = 0) = 1 



and for every n > 0, j G E and t > 0: 



(4) P{J n+l = j, X n+1 < t\(J k , X k ) k < n ) = P{J n+1 = j, X n+1 < t\(J n , X n )) = p Jnj F Jnj (t), 
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where p = (Pij)ijeE is a transition matrix and (Fij)ij eE is an array of distribution functions 



on M + = (0, +00). For more details on Markov Renewal processes see, for example, Limnios 
and Oprisan | ( pOOl ). We just recall that, under Assumptions ^ and Q: 



- the process ( J n ) n >o is a Markov chain, starting from jo, with transition matrix p, 

- the sojourn times (X n ) n > , conditionally on ( J n ) n >o, form a sequence of independent 
positive random variables, with distribution function Fj n _ 1 j n . 

We assume the functions are absolutely continuous with respect to the Lebesgue mea- 
sure with density fy. Hence, the likelihood function of data (j,x, ut) is: 

/r-l \ Hr>o) 

(5) L(j,x,ur) = lYlpjdi+Jjdi+i^i+i) I x^PjrkFjrkM, 

\i=0 J k&E 

where, for every x, F^ is the survival function: 

Fij(x) = 1 - Fij(x) = P(X n+1 > x\J n = i, J n+ i = j) . 

Furthermore, we assume that each inter-occurrence time density is a Weibull density 
with shape parameter and scale parameter 9 if 

/ \ an — 1 



(6) fij(?) = ir hr e y > x>o, ^>o, %>o. 

For conciseness, let ck = (ca^ij^E and = (0y)jj 6 £. 

In order to write the likelihood in a more convenient way, let us introduce the following 
natural statistics. We will say that the process visits the string if a visit to state % is 
followed by a visit to state j and denote: 

- x?j the time spent in state i at the p-th visit to string 

- Nij the number of visits to the string 

Then, assuming r > 0, from Equations ([5j and ^ we obtain the following expression for 
the likelihood function: 
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>< n 

i,keE 



(7) L(},x,v T \p,a,0)= 

i,k£E 

1 /N ik \ aife_1 f N ik V 

-(5»-{- (£)-}. 

Our purpose is now to perform a Bayesian analysis for p, a and 6 which allows us to 
introduce prior knowledge on the parameters. As we will show, this analysis is possible, 
using a Gibbs sampling approach, via the introduction of an auxiliary variable in order to 
treat easily the right-censored datum. 

4. Bayesian analysis 

4.1. The prior distribution. In our Bayesian analysis, we assume that p is independent 
on a. and 6. In particular, the rows of p are s independent vectors with Dirichlet distri- 
bution with parameters 71, • • • , ~f s and total mass ci, • • • c s , respectively. This means that, 
for i — 1, . . . , s, the prior density of the i-th row is 

Via 



7Tl;i(pa, • • • ,Pis) = ]r ,', J 



V 
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on T = {(pn, . . . ,p i8 )\ p i: j > 0, J2jPv = !} where 7» = (7a, • • • , lis) and a = £J =1 7^ . 

As far as a and is concerned, we will assume the 0j/s, given the a^-'s, are independent 
with generalized inverse Gamma density: 

In 1 N //i 1 \ a ijK( a ij) mij „-(l+m i7 -a i7 -) j hj{oiij) 1 „ 
7T2;y(%|a) = ^(fl^K) = r(m J) 6 ij X eX P j - : , J % > 0. 

where > and ^(a^) is a positive function of a^. In other terms, a priori ij - a,J , 
given ay, has a Gamma density with shape m,j and rate l/bij(aiij) > . In symbols 
0,j | «ij ~ GIG (rriij , &jj (otij ) , ay ) . 



Finally, we will assume the components of ex are independent and have a log-concave 



density n 3]i j with support bounded away from zero. As discussed in Gilks and Wild (1992), 



the log-concavity of 7^ j is necessary in the implementation of the Gibbs sampler (see also 



Berger and Sun (1993)), although adjustments exist for the non-log-concave case (see 
Gilks et al. (1995)). Conversely, a support bounded away from zero ensures the existence 
of the posterior moments of the 0y's. Anyway, we will specify TT^i j in detail in the next 
section. 

4.2. Elicitation of the hyperparameters. In this section we focus our attention on the 



prior of (9^,01^), for i,j fixed. Following an approach developed in Bousquet (2006), 
we provide a marginal prior distribution for the shapes along with a statistical 
justification of the choice of the generalized inverse Gamma prior. An interpretation of the 
hyperparameters is also provided. 

First of all, for the sake of clarity, let us drop the indices i,j in the notations and 
quantities introduced in the previous section. Suppose now that a vector of historical 
dataset x_ m = (x_ m , . . . , x_i) is available, i.e. a past "historical dataset" of m sojourn 
times in state i followed by a visit to state j. Moreover, let us accept to use as prior density 
for (a, 9) its posterior density given x_ m , when we start from the improper prior: 

OC0" 1 (l-^) C l (e > )l( a ; 



7r a, 



*>«o) 



with suitable c > and «o > 0. Then it is immediate to see that this action corresponds 
to assume as prior distribution for 9 given a the density: n 2 {9\a) = QIG(m, b(a, x_ TO ), a) 
and as prior density of a: 

with 6(a,x_ m ) = Y,?=i x -i and P[x~m) = "lEiM^-i)]" 1 - 

Remark 1. We notice that the prior for a, 9 we propose has a simple hierarchical structure: 

iT2(9\a) is a generalized inverse Gamma density with the first parameter equal to the size of 

an historical dataset from the same process, whereas the second parameter is a deterministic 
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function of a. Moreover, the improper prior (J8j) we start with depends on two parameters 
ao and c that can be chosen in a such way that the proper prior distribution of 9 derived 
from it has finite moments of sufficiently high order. Finally, it is easy to see that, if 
m > 1, then the prior of a is proper. 



4.3. Simplifying the prior. In the previous section we consider a prior with a simple 
hierarchical structure obtained as the posterior coming from a Bayesian inference on his- 
torical data. In order to simplify the expression of the posterior distributions and to ensure 
that posterior moments are finite, we replace the function b(a, x_ m ) = YlT=i x -ii * n ^(^Ict) 
and vr 3 (o;) by an easier convex function of a. More precisely, we replace such a function 
with b q (a) = t q [(l — q)~ x l m — where t q is a positive number and q G (0, 1). Then one 
can see that, for such a choice of b q (a), t q turns out to be the quantile of order q of the 
marginal distribution of the inter-occurrence times between a visit to state i followed by a 



visit to j. See Bousquet (2010). Hence, for a fixed q G (0, 1), we can estimate t q using the 



historical data. Denote by t q such an estimate and let 

b q (a^ m ) = i a q [(l-q)- l l m -l}- 1 . 
Thus, we propose for our Bayesian analysis the following priors: 

(10) 7r 2 (0|a) = giG(m, b g (a, x_ m ), a) 
and 

f ( - Y m - lnx_A 1 

(11) n 3 (a) oc a m 1 c (a - a ) c exp < -ma ( In t q - j \ l( Q ,> Q , o) . 

For m > 1, 7r 3 (a) is a proper prior if q is chosen in such a way that 

m 

and is always log-concave for any c > 0. Notice that such a choice is possible for almost 
all set of data. 
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Finally let us elicit the hyperparameters ao and c in order that the posterior second 
moment of 9 is finite. Since 

E(9 2 ) = E ( r( ^ ( "^ /a) [&,(«, *-m)] ^ < M(r(m - 2/a)) 

for a suitable constant K, then one can see that if ao = 2/m and c > 1, then E(# 2 ) < +00 
and hence too the posterior second moment of 9 is finite. As noticed above, ^(a) is a 
log-concave density for any c > 0. In particular, for c = m — 1, 7r 3 (a) is an ao-shifted 
Gamma density. 

4.4. Scarce prior information. The construction of the prior distribution of (a, 9) has 
to be modified for those pair of states between which no more than two transitions were 
observed, that is, m < 2. 

If m — 2 then qj = 1, which rules out decreasing hazard rates. In the absence of 
additional specific prior information, this is an arbitrary restriction, so an ao value smaller 
than 1 must be chosen. Then, the prior second moment of 9 is not finite anymore. For 
the posterior second moment to be finite we need a > 2/(2 + N), where N is the number 
of transitions between the two concerned states in the (current) sample. Thus the second 
moment of 9 can stay non-finite, even a posteriori, if 2/(2 + N) > ao- This would show 
that the data add little information for that specific transition. To avoid this, we may let 
«o = 2/3, corresponding to the smallest prior sample size such that ao < 1. The value of 
c can still be m — 1 = 1. 

If m = 1, the single historical observation x_i determines 6 g (a,x_ m ). As i q = x_i for 
any q, it seems reasonable to use q = 0.5, so x_i would represent the prior opinion on the 



median sojourn time. Since with m = 1 the argument of the exponent in 7T 3 (a) in (11) 

is zero so that 713(a) is improper, then we make it proper by restricting its support to an 

interval [accti]. The value a± = 10 is suitable for all practical purposes. As before, the 

choice ao = 2/m = 2 would be too much restrictive, so we select again ao = 2/3. The rule 

c = m — 1 would give c = in this case, so we keep c = 1. If m = 0, the prior information on 

the number of transitions is that there have been no transitions, but there is no information 
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on the sojourn time. Then we introduce a single random fictitious observation t q , which is 
uniformly distributed between the smallest and the largest sojourn times associated with 
all the transitions in the historical dataset, and we use t q to obtain b q (a,t q ). Thus we 
fall in the previous case by substituting m = 1 for m = and by assuming c = 1 and 
[a ,ai] = [2/3,10]. 

For clearness, we summarize the hyperparameter selection for priors (ph-dll} in Table 111 
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3 



Table 1. Hyperparameter selection as the prior sample size m varies 



5. Analysis of the central Northern Apennines sequence 

In this section we analyze the macroregion MR 3 sequence, using the semi-Markov model. 
Model fitting, model validation and an attempt at forecasting involves the following steps: 

(1) choose the historical dataset x_ m for the elicitation of the prior distribution; 

(2) model fit is assessed by comparing observed inter-occurrence times (grouped by 
transition) to posterior predictive intervals; 

(3) cross state-probabilities are estimated, as an indication to the most likely magnitude 
and time to the next event, given information up to the present time; 

(4) an interpretation in terms of slip predictable or time predictable model is provided. 

For the elicitation of the prior distribution we cut the observations into two parts, using 
the first part, called the historical data, for this purpose, and the second part, called the 
current data, for the likelihood. A look at the matrix of the observed transition frequencies 
in the whole dataset from year 1838 to 2002, displayed in Table |2^a), shows that our 195- 
event series contains enough data to allow the assignment of at least three events for any 
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transition to the historical dataset. This choice avoids the need to adjust the prior as 
summarized in Table [T] for m^- < 2. The cut-point corresponding to this strategy is the 
82nd event which occurred in 1916, with frequency matrix as in Table [2Fb). 



(a) (b) 
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65 


30 


17 
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14 
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32 


15 
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14 
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3 


15 


9 


4 


3 


4 


3 


3 



Table 2. Number of observed transitions in the whole (a) and historical (b) datasets. 



Let us consider the predictive check mentioned above. Fig. [3] shows posterior predictive 
95 percent probability intervals of inter-occurrence time for every transition, with the 
observed inter-occurrence times superimposed. These are empirical intervals computed 
by generating stochastic inter-occurrence times from their relevant distributions at every 
iteration of the Gibbs sampler. (The Gibbs sampler itself is completely illustrated in the 
appendix). Possible outliers, represented as triangles, are those times with Bayesian p- 
value (that is the predictive tail probability) less than 2.5 percent. The percentage of these 
extreme points, along with the predictive expected value of the inter-occurrence times, 
is showed in Table |3j Most points fall within the predictive intervals. Those few inter- 
occurrence times which are really extreme, such as the small values observed in transitions 
(1, 1), (2, 1), (1, 2) and the large one in transition (1, 3) match unsurprisingly the outlying 
points in the corresponding qq-plots in Fig. [2j (In some qq-plots there are more outlying 
points than in Fig. [3] because the qq-plots comprise the entire series). This fact could be 
regarded as a lack of fit of the Weibull model, but we are more inclined to attribute it to 
an imperfect assignment of some events to the macroregion or to an insufficient filtering of 
secondary events. 

The shape parameters are particularly important as they reflect an increasing hazard 
if larger than 1, a decreasing hazard if smaller than 1 and a constant hazard if equal to 1. 
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FIGURE 3. Posterior predictive 95 percent credible intervals of the mean inter- 
occurrence times in days with actual times denoted by (blue) solid dots. Suspect 
outliers are denoted by (red)-pointing triangles. The (green) dotted line shows the 
posterior mean and the (violet) dashed line the posterior median. 



1 2 3 

1 253.28 (10.5%) 272.27 (18.8%) 340.98 (7.7%) 

2 273.45 (16.7%) 314.67 (14.3%) 293.35 (3.3%) 

3 204.48 ( 0.0%) 237.54 (16.7%) 312.73 (0.0%) 

Table 3. Summaries of the posterior distributions of the mean inter- 
occurrence times for each transition. Posterior means are followed by (per- 
centage of points having posterior p- value less than 2.5 percent). 
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1 


2 


3 


1 


1.182 (0.102) 


1.981 (0.252) 


0.828 (0.103) 


2 


1.318 (0.173) 


1.218 (0.232) 


1.716 (0.473) 


3 


1.026 (0.179) 


0.996 (0.150) 


2.098 (0.618) 



Table 4. Summaries of the posterior distributions of the shape parameter 
a. Posterior means are followed by (standard deviations). 



Table [4] displays the posterior means of these parameters (along with their posterior stan- 
dard deviations). Finally Table [5] shows the posterior means of the transition probabilities. 





1 


2 


3 


1 


0.57 


0.27 


0.16 


2 


0.58 


0.28 


0.14 


3 


0.52 


0.32 


0.16 



Table 5. Posterior means of the transition matrix p. 



Cross state-probability plots are an attempt at predicting what type of event and when 
it is most likely to occur. A cross state-probability (CSP) P^ Ax represents the probability 
that the next event will be in state j within a time interval Ax under the assumption that 
the previous event was in state % and to time units have passed since its occurrence: 



( 12 ) P Z\^ = P ( J n+l = j, X n+1 <t + Ax\J n = i, X n+1 > t 







Pij (Fijito) -F ij (t + Ax)) 



Y.hdEPihFih{t ) 

Fig. [4] displays the CSPs with time origin on 31 December 2002, the closing date of the 



CPTI04 (2004) catalogue. At this time, the last recorded event had been in class 2 and 
had occurred 965 days earlier (so to is about 32 months). From these plots we can read 
out the probability that an event of any given type will occur before a certain number of 
months. For example, after 24 months, the sum of the mean CSPs in the three graphs 
indicates that the probability that an event will have occurred is more than 90 percent, 
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with a larger probability assigned to an event of type 2, followed by type 1 and type 3. 
The posterior means of the CSPs are also reported in Table |6j 



90% posterior predictive intervals of CSP from 2 to 1 
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90% posterior predictive intervals of CSP from 2 to 2 
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90% posterior predictive intervals of CSP from 2 to 3 
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FIGURE 4. Posterior mean and median of CSPs with time origin on 31 December 
2002 up to 48 months ahead, along with 90 percent posterior credible intervals. 
Transitions are from state 2 to state 1, 2 or 3 (first to third panel, respectively). 





1 Month 2 Months 


3 Months 


4 Months 


5 Months 6 Months 


Year 2 Years 


3 Years 


4 Years 


to 1 
to 2 
to 3 


0.070 0.122 
0.059 0.105 
0.014 0.024 


0.171 
0.150 
0.034 


0.210 
0.187 
0.041 


0.243 0.270 
0.221 0.249 
0.048 0.054 


0.364 0.410 
0.363 0.444 
0.075 0.088 


0.418 
0.468 
0.092 


0.420 
0.477 
0.093 




Table 6. CSPs 


with time 


origin on 


31 December 2002, 


as represented 


in Fig. 
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The predictive capability of our model can be assessed by marking the time of the next 
event on the relevant CSP plot. In our specific case, the first event in 2003 which can 
be assigned to the macroregion MR 3 happened in the Forli area on 26 January and was 
of type 1, with a CSP of 7 percent. This is a low probability, but a single case is not 
enough to judge our model, which would be a bad one if repeated comparisons did not 
reflect the pattern represented by the CSPs. Therefore we repeated the same comparison 
by re-estimating the model using only the data up to 31 December 2001, 31 December 
2000, and so on backwards down to 1992. For this particular exercise the historical data 
for prior elicitation are such that the frequency of each transition is at least two (not three 
as done so far) and the hyperparameters are those on the second column of Table [I] The 
reason is that, by going backwards, the only (3, 3) transition is removed from the current 
data set, so that the CSP for (3, 3) would rely on prior information only, resulting in a too 
large posterior variance of the parameters. The results are shown in Table [7| The boxed 
numbers correspond to the observed events and it is a good sign that they do not always 
correspond to very high or very low CSPs, as this would indicate that events occur too 
late or too early compared to the estimated model. 

An exception is represented by the last four arrays in the table, which are all related to 
the period from the last event in July 1987 to the event in October 1996, during which 
no events with M w > 4.5 occurred. This is the longest inter-occurrence time in the whole 
series and is the outlying point in the top right panel of Fig. [2] and the associated transition 
is a relatively rare one, so it is not surprising that the boxed CSPs are small. 

The examination of the posterior distributions of transition probabilities and of the 
predictive distributions of the inter-occurrence times can give some insight into the type of 
energy release and accumulation mechanism. We consider two such mechanisms, the time 
predictable model (TPM) and the slip predictable model (SPM). 

In the TPM, when a maximal energy threshold is reached, some fraction of it (not always 
the same) is released; therefore, the waiting time for the next event is longer if the current 
event is stronger. So, the waiting time distribution depends on the current event type, but 
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end of catalogue: 31/12/2001; previous event type:2; waiting time: 600 days 





1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


Year 


392 days 


2 Years 


3 Years 


4 Years 


to 1 


0.053 


0.095 0.135 


0.169 0.199 0.224 


0.318 


0.373 


0.375 


0.384 


0.387 


to 2 


0.034 


0.063 0.092 


0.118 0.143 0.166 


0.272 


0.383 


0.388 


0.432 


0.456 


to 3 


0.017 


0.031 0.043 


0.054 0.063 0.070 


0.098 


0.115 


0.115 


0.119 


0.121 






end of catalogue: 31/12/2000; previous event type: 2; 


waiting time: 235 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


Year 


2 Years 


757 days 


3 Years 


4 Years 


to 1 


0.061 


0.112 0.162 


0.205 0.245 0.280 


0.417 


0.502 


| 0.505 


0.518 


0.522 


to 2 


0.022 


0.040 0.059 


0.076 0.092 0.106 


0.176 


0.248 


0.252 


0.279 


0.293 


to 3 


0.017 


0.031 0.047 


0.061 0.074 0.085 


0.132 


0.159 


0.159 


0.164 


0.165 






end of catalogue: 31/12/1999; previous event type: 1; 


waiting time: 177 days 






1 Month 


2 Months 3 Months 


4 Months 130 days 5 Months 


6 Months 


Year 


2 Years 


3 Years 


4 Years 


to 1 


0.062 


0.115 0.166 


0.210 0.222 0.251 


0.286 


0.430 


0.525 


0.544 


0.548 


to 2 


0.038 


0.074 0.110 


0.142 1 0.151 1 0.172 


0.197 


0.283 


0.305 


0.305 


0.305 


to 3 


0.013 


0.023 0.033 


0.042 0.045 0.050 


0.058 


0.091 


0.122 


0.135 


0.140 






end of catalogue: 31/12/1998; previous event type: 3; 


waiting time: 280 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


188 days 


Year 


2 Years 


3 Years 


4 Years 


to 1 


0.062 


0.109 0.153 


0.189 0.221 0.247 


0.252 


0.339 


0.389 


0.400 


0.402 


to 2 


0.037 


0.067 0.095 


0.119 0.141 0.160 


0.164 


0.234 


0.290 


0.307 


0.313 


to 3 


0.037 


0.068 0.100 


0.126 0.150 0.170 


0.175 


0.238 


0.268 


0.274 


0.276 






end of catalogue: 31/12/1997; previous event type: 3; 


waiting time: 442 days 






1 Month 


2 Months 85 days 


3 Months 4 Months 5 Months 


6 Months 


Year 


2 Years 


3 Years 


4 Years 


to 1 


0.074 


0.131 0.176 


0.184 0.227 0.264 


0.294 


0.402 


0.463 


0.475 


0.479 


to 2 


0.054 


0.096 0.131 


0.138 0.173 0.205 


0.232 


0.343 


0.428 


0.455 


0.466 


to 3 


0.010 


0.016 0.021 


0.021 0.025 0.028 


0.030 


0.037 


0.041 


0.042 


0.043 






end of catalogue: 31/12/1996; previous event type: 3; 


waiting time: 77 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


Year 


450 days 


2 Years 


3 Years 


4 Years 


to 1 


0.073 


0.131 0.186 


0.231 0.271 0.304 


0.422 


0.447 


0.483 


0.493 


0.496 


to 2 


0.043 


0.075 0.106 


0.132 0.155 0.174 


0.248 


0.268 


0.300 


0.314 


0.319 


to 3 


0.012 


0.027 0.049 


0.076 0.105 0.129 


0.172 


0.175 


0.178 


0.179 


0.179 






end of catalogue: 


31/12/1995; previous event type: 1; 


waiting time: 3100 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


288 days 


Year 


2 Years 


3 Years 


4 Years 


to 1 


0.173 


0.304 0.420 


0.511 0.589 0.651 


0.798 


0.861 


0.964 


0.982 


0.985 


to 2 


0.001 


0.002 0.002 


0.003 0.003 0.004 


0.004 


0.004 


0.005 


0.005 


0.005 


to 3 


0.001 


0.003 0.004 


0.004 0.005 0.005 


0.007 


0.007 


0.008 


0.008 


0.008 






end of catalogue: 


31/12/1994; previous event type: 1; 


waiting time: 2735 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


Year 


653 days 


2 Years 


3 Years 


4 Years 


to 1 


0.169 


0.295 0.410 


0.502 0.580 0.643 


0.857 


0.954 


0.964 


0.982 


0.986 


to 2 


0.001 


0.002 0.002 


0.003 0.003 0.004 


0.005 


0.005 


0.005 


0.005 


0.005 


to 3 


0.001 


0.002 0.003 


0.004 0.004 0.005 


0.007 


0.007 


0.008 


0.008 


0.008 






end of catalogue: 


31/12/1993; previous event type: 1; 


waiting time: 2370 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


Year 


2 Years 


1018 days 


3 Years 


4 Years 


to 1 


0.166 


0.289 0.403 


0.495 0.573 0.636 


0.854 


0.964 


0.981 


0.982 


0.986 


to 2 


0.001 


0.002 0.003 


0.003 0.004 0.004 


0.005 


0.005 


0.005 


0.005 


0.005 


to 3 


0.001 


0.002 0.003 


0.003 0.004 0.004 


0.006 


0.007 


1 0.007 


0.007 


0.007 






end of catalogue: 


31/12/1992; previous event type: 1; 


waiting time: 2005 days 






1 Month 


2 Months 3 Months 


4 Months 5 Months 6 Months 


Year 


2 Years 


3 Years 


1383 days 


4 Years 


to 1 


0.161 


0.283 0.395 


0.486 0.564 0.628 


0.849 


0.963 


0.983 


0.986 


0.987 


to 2 


0.001 


0.002 0.003 


0.003 0.004 0.004 


0.005 


0.006 


0.006 


0.006 


0.006 


to 3 


0.001 


0.002 0.002 


0.003 0.003 0.004 


0.005 


0.006 


0.006 


0.006 


0.006 



TABLE 7. CSPs as the end of the catalogue shifts back by one- year steps. The 



numbers in boxes are the probability that the next observed event has occurred 
at or before the time when it occurred and is of the type that has been observed. 
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not on the next event type, that is, we expect F^t) = Fi.(t), j = 1, 2, 3. The strength of an 
event, does not depend on the strength of the previous one, because every time the same 
energy level has to be reached for the event to occur, so we expect pij = p.j, j = 1, 2, 3, that 
is, a transition matrix with equal rows (which is indeed the case, as seen from Table |5]). 



The CSPs ( 12 ) would factorize as follows, 

( , pij = Pa {Fijjto) ~ Fjjjto + Ax)) = p., (FUM - Fi.(t + Ax)) 

1 ' tolAx J2 h e E PihF ih (t ) F,(t ) 

so that, given i, they are proportional to each other as j = 1,2,3 for any Ax. 

In the SPM, after an event energy falls to a minimal threshold and increases until the 
next event, where it starts to increase again from the same threshold. Here, the magnitude 
of an event depends on the length of the waiting time, but not on the magnitude of the 
previous one, because energy always accumulates from the same threshold. In this case 
again p^ = p.j, but Fij(t) = F.j(t), so 



_ /,,,• (/•:,•(/„) - /•;,•(/„ • A,-)) 



If this is the case CSPs, given j would be equal to each other as i — 1, 2, 3 for any Ax. 
As the a posteriori transition probability matrix has essentially equal rows, the CSPs 



should be examined with regard to the properties ( 13 ) and ( 14 ) to decide whether our se- 
quence is closer to a TPM or to an SPM. An additional feature that can help discriminate 
is the tail of the waiting time distribution: for a TPM, the tail of the waiting time distribu- 
tion is thinner after a weak earthquake than after a strong one; for an SPM, the tail of the 
waiting time is thinner before a weak earthquake than before a strong one. However this 
additional job proves unnecessary, after examination of Figures [5] and [6j which represent 
the posterior means of the ratios of the CSPs as a function of Ax for to = 0. In Fig. |5j the 
ratio of two CSPs -say -P t ^i Aa ./ Pl^ x ~ starts to be equal to Pij/pik only when Ax is large 
enough for the survival functions Fij(to + Ax) and Fik(to + Ax) to be close to zero. A 
similar behaviour is observed in Fig. [6j The conclusion is that neither a SPM nor a TPM 
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is supported by the posterior distributions of the parameters, because of the behavior of 
the waiting times. 



6. Concluding Remarks 

Whe have presented a complete methodology for the Bayesian inference of a semi-Markov 
process, from the elicitation of the prior distribution, to the computation of posterior sum- 
maries, including a guidance for its JAGS implementation. In particular, we have examined 
in detail the elicitation of the joint prior density of the shape and scale parameters of the 
Weibull-distributed holding times (conditional on the transition between two given states), 
deriving in a natural way a specific class of priors, along with a method for the determi- 
nation of hyperparameters based on "historical data" and moment existence conditions. 
This framework has been applied to the analysis of seismic data, but it can be adopted for 
inference on any system for which a Markov renewal process model is plausible. A possible 
and not- yet explored application is the modelling of voltage sags (or voltage dips) in power 
engineering: the state space would be formed by different classes of voltage, starting from 
voltage around its nominal value, down to progressively deeper sags. In the engineering 
literature, the dynamic aspect of this problem is in fact disregarded, while it could help 
bringing additional insight into this phenomenon. 

With regard to the seismic data analysis, other uses of our model can be envisaged. The 
model can be applied to areas with a less complex tectonics, such as Turkey, by replicating 
for example Alvarez's anaysis. Outliers, such as those appearing in Figure [3j could point 
out to events whose assignment to a specific seismogenetic source should be re-discussed. 
The analysis of earthquake occurrence can support decision making related to the risk of 
future events. We have not examined this issue here, but a methodology is outlined by 



Cano et al. (2011) 



A final note concerns the more recent Italian seismic catalogue CTPI11 (2011 ), including 



events up to the end of 2006. Every new release of the catalogue involves numerous changes 
in the parameterization of earthquakes; as the DISS event classification by macroregion is 
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Ratios of couples of CSPs from 1 



Ratios of couples of CSPs from 2 



-«- ratio of CSP (1,1) to (1,2) 
• ratio of CSP (1,1) to (1,3) 
-v - ratio of CSP (1,2) to (1,3) 



_o — o O O O O- 



i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 

1M 3M 5M 1Y 2Y 3Y 4Y 5Y 6Y 7Y 8Y 9Y 



o- ratio of CSP (2,1) to (2,2) 
• ratio of CSP (2,1) to (2,3) 
-v - ratio of CSP (2,2) to (2,3) 



I — 
1M 



~i 1 1 1 1 1 1 1 1 1 r 

1Y 2Y 3Y 4Y 5Y 6Y 7Y 8Y 9Y 



(a) 



(b) 



Ratios of couples of CSPs from 3 



— i — 

1M 



-o- ratio of CSP (3,1) to (3,2) 
• ratio of CSP (3,1) to (3,3) 
-V - ratio of CSP (3,2) to (3,3) 



~i 1 1 1 1 1 1 1 1 1 r 

1Y 2Y 3Y 4Y 5Y 6Y 7Y 8Y 9Y 



(c) 

FIGURE 5. Posterior means of ratios of CSPs ■Fgj'^ a ./-Foj^\ a . with time origin on 
up to 10 years ahead. Transitions are from state 1, 2 and 3 to state 1, 2 or 3. 
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Ratios of couples of CSPs to 1 



Ratio of couples of CSPs to 2 



o- ratio of CSP (1,1) to (2,1) 
• ratio of CSP (1,1) to (3,1) 
-V - ratio of CSP (2,1) to (3,1) 



w ¥¥¥¥¥¥¥¥ 

■ 

_V Q o O O O O O O 0- 



i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 

1M 3M 5M 1Y 2Y 3Y 4Y 5Y 6Y 7Y 8Y 9Y 





o 

o 


-«- ratio of CSP (1,2) to (2,2) 
• ratio of CSP (1,2) to (3,2) 
-v - ratio of CSP (2,2) to (3,2) 


• 
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V V V V V V V 


° V 
V 

V 

v o 

V 

V 

• 




• 

o 

• 




° • 




• 




• 





1Y 2Y 3Y 4Y 5Y 6Y 7Y 8Y 9Y 



(a) 



(b) 



Ratio of couples of CSPs to 3 



o- ratio of CSP (1 ,3) to (2,3) 
• ratio of CSP (1 ,3) to (3,3) 
-V - ratio of CSP (2,3) to (3,3) 



— I — 
1M 



S — S — S — fi — sz — a — v — sz — sz_ 

~i — i — i — i — i — i — i — i — i — i — r 

1Y 2Y 3Y 4Y 5Y 6Y 7Y 8Y 9Y 



(c) 



FIGURE 6. Posterior means of ratios of CSPs Pq\^ x / Pq\& x with time origin on 

up to 10 years ahead. Transitions are from state 1, 2 and 3 to state 1, 2 or 3. 
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not yet available for events in this catalogue we could not use this more recent source of 
data. 



Appendix A. Gibbs sampling 

Here we derive the full conditional distributions needed to implement a Gibbs sampling 
and give indications on its JAGS (Plummer 2010) implementation. 

A.l. Full conditional distributions. Without loss of generality we can assume that 
the last holding time is censored i.e. ut > 0. In this case in order to obtain simple full 
conditional distributions and then an efficient Gibbs sampling, we introduce the auxiliary 
variable j T+ i which represents the unobserved state following j T . Thus the state space of 
the Gibbs sampler is (p,a,6,j T+ i) and the following full likelihood derived from (m): 



(15) L(j,x,u r ,j T+ i|p,O!,0) = Yl Pfk ikx 



X 



n 

i.keE 



.Ma 



■)<H k N ik ( II X ik ) 

7 ik \p=i J 



«ifc-l 



j N ik 

' f) a ik 
ik „_i 



x exp 

p=l 

x (ftVi T+ iexp 



x 



P \Ctik 

ik) 



x 3t3 t + 1 



is multiplied by the prior and used to determine the full conditionals. For every i and j let 



PH) 

a(-ii) 



e 



M ij (a ij ) 



the transition matrix p without z-th row, 
(a hk , h,k e E, (h,k) ^ , 
{6m,, h,keE, (h,k) ^ , 



Hi) 



N i:j = Nij + l(ov,i T+1 )=(ij)), Ni=[Nij, j = 1, . 



EW aw + ^i((jV,j. + i) = (^,j)), 
p=l 



^ = n 

p=l 



4 
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Furthermore, let rriij be the number of visits to the string in the historical dataset, 
x~ p , for p = l,...TOij, denote the corresponding inter-occurrence times and x_ mij . = 

/ —rriij _]\ 
V-^ij j • • • j ^ij ) ■ 

The following result on the full conditional distributions of the Gibbs sampling holds. 



Proposition A.l. Let the prior on (p, a, 6) be as in Section^ Then 

(a) the conditional distribution ofpi, given}, x, ut, jV+i,P(-«), & andO, is a Dirichlet 
distribution with parameter Nj + 7, ; 

(b) the conditional distribution of 9^ , given}, x, ut, j T +i, P, ot and Q(-ij) is an 
inverse- Gamma distribution of parameters rriij + and b q (aij, x_ my ) + M^cty) ; 

(c) the conditional density of a^, given }, x, ut, jV+i, P, and is log-concave 
and, up to a multiplicative constant, is equal to: 

(i6) x ( Cii ti er r -»>f> exp / jiteaH^ 

where t l j.. is determined from Table 

((i) t/ie conditional density of the unseen state J T +i, given}, x, «t, p, a and ; zs 



PjVj ex P 



E, 6S ^ex P {-(^-) Q3TJ }' 

(e) i/m^ > thent^ is a known constant (see Table 1). If rriij = 0, then the conditional 
distribution of t l J.. given}, x, ut, p, cx and 6 is a double-truncated Weibull with shape 
aij + 1, scale 8ij and the truncation points given by the smallest and the largest inter- 
occurrence times included in all the historical dataset. 



Proof. Notice that the row p^, given the data and j' T +i is conditionally independent on 
(p(-i), oc, 0). Hence 



7r(pi|j,x,«T, j T+ i,p^i),a,e) oc L(j,x,w T , j T+1 \p,a,0) x 7ri ;i (pi) oc x JJpJ 

j&E j&E 

24 



and point (a) of the thesis follows. 

With regard to the full conditional distribution of % given (j, x, ut, jV+i, P, a , we 
have 



7T 2 : 



|j,X,U T , J r +l,P,Q!,0(-y)) OcL(j,X,M T , j T+1 |p,Q!,0) X 7T 2 



ij{uij\0>ijj 



II 4 U >< II 



i,k£E 



x exp 



ttife-1 



M ik( a ik) 

f, a ik 

x e « 



x 



Tfrn 



^5 (^ij 5 ~X-—m,ij , 
'J 



oc exp 



6 q (aij,x_ m<J .) + A/ij^ajj-; 



As one can see, the last function is the kernel of an inverse Gamma distribution with 

parameters + % and b q (aij, x_ mij .) + M;j(ajj) and point (6) of the thesis follows. 

A similar reasoning yelds a posterior distribution 7r 3;i j(ajj|j, x, ut, jV+i, P, Qt(-y), 0) which 



is proportional to (16). Furthermore, concerning its log-concavity, from (16) it follows that 



log7r(a y |j, x, u T , j T +i, P, 0) = 

= log ^(atij) + otij log(Cyt^) - (rriij + Ny) log 6^- 



+ {Nij + 1) log 



a; 



All terms of this expression are log-concave. In particular ir^-ij is log-concave by con- 
struction and the last term is a sum of concave functions of kind g(z) t— > —z aij . Hence 
the log-concavity follows from the property that the sum of concave functions is concave. 
Finally (El) implies the following: 

P{Jt+i = j, X T+1 > u T \ j, x, p, cx, 0) 



P(Jt+i = j\ j,x,u T ,p,oi,0) 



Y.keE P ( J T+i = K x r+i > u T \ j,x,p,cx,0) 

Pjrj eX P 



9 iTj 
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ut 



a irk 



The same type of argument can be used for proving (e). □ 



A. 2. JAGS implementation. Proposition A.l implies that JAGS should be able to run 



an exact Gibbs sampler. The model description we have adopted in the JAGS language 



based on the full likelihood (15). It is not important that the model description matches 
the actual model which generated the data as long as the full conditional distributions 
remain unchanged. 

In particular, for any i, j T and Ylk -^ik can be considered as fixed. Then, given pj, 
(Nn,...,N is ) is regarded as a sample from a multinomial distribution with probability 
vector pi and J2k N%k trials; furthermore, the factor Pj T j T+1 in the full likelihood is regarded 
as the prior for J T +\ as it ranges from 1 to s; finally, upon multiplication by the Dirichlet 



prior for pi we re-obtain the correct full conditional distribution of Proposition A.l (a). 

The uncensored inter-occurrence times are described as Weibull distributed with shape 
and rate parameters that depend on j. The right-censored observation ut is handled by 
a special instruction provided by the JAGS language, which should represent the actual 
likelihood term correctly, so that the correct form of the full conditional of J T +\ is preserved. 

The description of the prior distributions is routine. We have mentioned the Dirichlet 
prior for pj; ctij has a shifted Gamma distribution if > 1, as specified at the end of 
Section 



4.3 



:= l/9™j j is Gamma distributed, conditionally on a^, for any value of my 
and t^ ; i 1 ^. is uniformly distributed if = 0, otherwise it is constant. The only nontrivial 
case arises for the prior of if < 1, which reduces to ^3(0^) oc (1 — a /a) l ao < a < ai 



see Section 4.4). This distribution is coded using the so-called zeros trick: a fictitious 



zero observation from a Poisson distribution with mean = —]nTr(ctij) is introduced; a 
uniform prior over [ao,ai] is assigned to a^; then the likelihood factor corresponding to 
the zero observation is exp(— 0), which multiplied by the uniform prior gives the desired 
factor in the joint distribution of the data and the parameters. 
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